#include "vainshtein.h"

#define SIGN(a) (a>0)?1:0
#define EM_LS  GSL_CONST_MKSA_SPEED_OF_LIGHT // 2.99792458e8	//скорость света
#define M_MAX 80 // number of gn for fi&xi

Vainshtein::Vainshtein()
{
    for(int j=0;j<9; j++) RES[j]=0;
}

Vainshtein::Vainshtein(Integrals *ft)
{
    waveType=ft->waveType;
    mI=ft->mI; nI=ft->nI;
    a=ft->a;    k=ft->k;
    R=ft->R; teta=ft->teta; fi=ft->fi; Z=ft->Z; Dist=ft->Dist;
    Shield=a*ft->difQ[3];
    for (int j=0; j<=3; j++) difQ[j]=ft->difQ[j];
    for (int j=1; j<=M_MAX; j++)
    {
        VMN[j-1]=gsl_sf_bessel_zero_Jnu(mI,j);
        MMN[j-1]=zerosJs(mI,j);
        if (Shield)
        {
            EMN[j-1]=zerosCoaxE(mI,j);
            HMN[j-1]=zerosCoaxH(mI,j);
        }
    }
    for (int j=0; j<40; j++) { Rln[j]=0; RlnQ[j]=0; Tln[j]=0; TlnQ[j]=0; }
    for(int j=0;j<9; j++) RES[j]=0;
    Kappa=k*a;
    GSL_SET_COMPLEX(&eps,1,0);
    if (VMN[nI-1] < Kappa) Gamma=sqrt(pow(Kappa,2)-pow(VMN[nI-1],2)); //: Gamma=1I*sqrt(pow(VMN[nI-1],2)-pow(Kappa,2));
    if (MMN[nI-1] < Kappa) GammaT=sqrt(pow(Kappa,2)-pow(MMN[nI-1],2)); //: GammaT=1I*sqrt(pow(MMN[nI-1],2)-pow(Kappa,2));
    Delta=compSqrt(gnMultiply(Kappa)/gnMultiplyT(Kappa))*eeDelta(Kappa);;//*eSw(Kappa)/eSwT(Kappa);
}

void Vainshtein::coeffs(QProgressBar *progress_int)
{
    progress_int->setMaximum(M_MAX);
    progress_int->setValue(0);	// progress-bar
    double dV,dVT,wn,wnt,gn,gnt,dW,dWT;   //diskrets for derivative
    __complex__ double fiPlu,xiPlu;
    __complex__ double fiPluS=0,fiMinuS=0,xiPluS,xiMinuS;
    gsl_complex expIndex,gslBuf;
    __complex__ double SUM[4];
    for(int j=0;j<9; j++) RES[j]=0;	//Инициализация
    switch (waveType)
    {
    case 1 :    // E-type
        fiPlu=(-1)*sqrt((Kappa+Gamma)*Gamma/VMN[nI-1])*compSqrt(hankel(mI,VMN[nI-1])*besselS(mI,VMN[nI-1])*gnMultiply(Gamma))*eSw(nI-1); // for frac
        numCoef[1]=0;
        while(VMN[numCoef[1]] < Kappa) numCoef[1]++;
        for (int j=1; j<=numCoef[1]; j++)
        {
            progress_int->setValue(j);	// progress-bar
            gn=sqrt(pow(Kappa,2)-pow(VMN[j-1],2));
            if (MMN[j-1] < Kappa) gnt=sqrt(pow(Kappa,2)-pow(MMN[j-1],2));
            wn=gn/a;
            wnt=gnt/a;
            //coefficients
            Rln[j-1]=fiPlu/sqrt((Kappa-gn)*gn/VMN[j-1])*compSqrt(gnMultiply(gn)/(hankel(mI,VMN[j-1])*besselS(mI,VMN[j-1])))*eSw(j-1);
            if (MMN[j-1] < Kappa) Tln[j-1]=fiPlu/sqrt((Kappa+gnt)*gnt/MMN[j-1])*compSqrt(gnMultiplyT(gnt)/(hankelS(mI,MMN[j-1])*besselS2(mI,MMN[j-1])))*eSwT(j-1);
            Rln[j-1]*=getRln(gn);
            Tln[j-1]*=getTln(gnt);
            GSL_SET_COMPLEX(&expIndex,__real__(1I*wn*Z),__imag__(1I*wn*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[0]=Rln[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(1I*wnt*Z),__imag__(1I*wnt*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[1]=Tln[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(-1*1I*wn*Z),__imag__(-1*1I*wn*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[2]=RlnQ[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(-1*1I*wnt*Z),__imag__(-1*1I*wnt*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[3]=TlnQ[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            //add to result
            RES[3]=SUM[1]+SUM[3];  //i_Z
            RES[6]=SUM[0]+1I*mI*SUM[1]+SUM[2]-1I*mI*SUM[3];  //i_fi
        }
        break;
    case 2 :    // H-type
        xiPlu=(-1)*sqrt((Kappa-GammaT)*GammaT/MMN[nI-1])*compSqrt(hankelS(mI,MMN[nI-1])*besselS2(mI,MMN[nI-1])*gnMultiplyT(GammaT))*eSwT(nI-1); // for frac
        numCoef[2]=0;
        while(MMN[numCoef[2]] < Kappa) numCoef[2]++;
        for (int j=1; j<=numCoef[2]; j++)
        {
            progress_int->setValue(j);	// progress-bar
            if (VMN[j-1] < Kappa) gn=sqrt(pow(Kappa,2)-pow(VMN[j-1],2)); else gn=0;
            gnt=sqrt(pow(Kappa,2)-pow(MMN[j-1],2));
            wn=gn/a;
            wnt=gnt/a;
            //! coefficients
            Rln[j-1]=xiPlu/sqrt((Kappa+gnt)*gnt/MMN[j-1])*compSqrt(gnMultiplyT(gnt)/(hankelS(mI,MMN[j-1])*besselS2(mI,MMN[j-1])))*eSwT(j-1);
            //Rln[j-1]=eSwT(j-1);
            if (VMN[j-1] < Kappa) Tln[j-1]=xiPlu/sqrt((Kappa-gn)*gn/VMN[j-1])*compSqrt(gnMultiply(gn)/(hankel(mI,VMN[j-1])*besselS(mI,VMN[j-1])))*eSw(j-1);
            //Tln[j-1]=eSw(gn);
            Rln[j-1]*=getRlnT(gnt);
            Tln[j-1]*=getTlnT(gn);
            //! sums with exponents
            GSL_SET_COMPLEX(&expIndex,__real__(1I*wnt*Z),__imag__(1I*wnt*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[0]=Rln[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(1I*wn*Z),__imag__(1I*wn*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[1]=Tln[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(-1*1I*wnt*Z),__imag__(-1*1I*wnt*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[2]=RlnQ[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            GSL_SET_COMPLEX(&expIndex,__real__(-1*1I*wn*Z),__imag__(-1*1I*wn*Z));
            gslBuf=gsl_complex_exp(expIndex);
            SUM[3]=TlnQ[j-1]*(GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf));
            //add to result
            RES[3]=SUM[1]+SUM[3];  //i_fi
            RES[6]=1I*mI*SUM[0]+SUM[1]-1I*mI*SUM[2]+SUM[3];  //i_Z
        }
        break;
    }
//    RES[3]*=B*sin(mI*fi);
//    RES[6]*=B*cos(mI*fi);
}

double Integrals::besselS(int m, double x)
{
    return m/x*gsl_sf_bessel_Jn(m,x)-gsl_sf_bessel_Jn(m+1,x);
}

double Integrals::besselS2(int m, double x)
{
    return (pow(m/x,2)-1)*gsl_sf_bessel_Jn(m,x)-besselS(m,x)/x;
}

__complex__ double Vainshtein::hankel(int m, double x)
{
    return gsl_sf_bessel_Jn(m,x)+1I*gsl_sf_bessel_Yn(m,x);
}

__complex__ double Vainshtein::hankelS(int m, double x)
{
    return besselS(m,x)+1I*(m/x*gsl_sf_bessel_Yn(m,x)-gsl_sf_bessel_Yn(m+1,x));
}

__complex__ double Vainshtein::compSqrt(__complex__ double val)
{
    gsl_complex gslBuf;
    GSL_SET_COMPLEX(&gslBuf,__real__(val),__imag__(val));
    gslBuf=gsl_complex_sqrt(gslBuf);
    val=GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf);
    return val;
}

__complex__ double Vainshtein::getRln(double gN)
{
    __complex__ double Rln;
    if ((Gamma+gN)==0) return 1;
    Rln=sqrt((Kappa+Gamma)/(Kappa-gN))/(Gamma+gN);
    Rln*=1-1/(1+1/Delta/Delta)*(Gamma+gN)/(Kappa+Gamma)/(Kappa+gN);
    return Rln;
}

__complex__ double Vainshtein::getTln(double gNt)
{
    __complex__ double Tln;
    Tln=Delta*(-2)/(1+Delta*Delta)/sqrt((Kappa+Gamma)*(Kappa-gNt));
    return Tln;
}

__complex__ double Vainshtein::getRlnT(double gNt)
{
    __complex__ double RlnT;
    if ((GammaT+gNt)==0) return 1;
    RlnT=sqrt((Kappa-gNt)/(Kappa+GammaT))/(GammaT+gNt);
    RlnT*=Delta*Delta/(1+Delta*Delta)*2*Kappa*(GammaT+gNt)/(Kappa-GammaT)/(Kappa-gNt)+1;
    return RlnT;
}

__complex__ double Vainshtein::getTlnT(double gN)
{
    __complex__ double TlnT;
    TlnT=pow(Kappa,2)/sqrt((Kappa+GammaT)*(Kappa-gN));
    TlnT*=-2*Delta/(1+Delta*Delta)/(Kappa-GammaT)/(Kappa+gN);
    return TlnT;
}

__complex__ double Vainshtein::gnMultiply(__complex__ double ww)
{
    __complex__ double Mul,gn; Mul=1+0I;
    for (int j=0; j<M_MAX; j++)
    {
        (VMN[j] < Kappa) ? gn=sqrt(pow(Kappa,2)-pow(VMN[j],2)) : gn=1I*sqrt(pow(VMN[j],2)-pow(Kappa,2));
        if (ww==gn) Mul*=(gn+ww); else Mul*=(gn+ww)/(gn-ww);
    }
    return Mul;
}

__complex__ double Vainshtein::gnMultiplyT(__complex__ double ww)
{
    __complex__ double Mul,gnt; Mul=1+0I;
    for (int j=0; j<M_MAX; j++)
    {
        (MMN[j] < Kappa) ? gnt=sqrt(pow(Kappa,2)-pow(MMN[j],2)) : gnt=1I*sqrt(pow(MMN[j],2)-pow(Kappa,2));
        if (ww==gnt) Mul*=(gnt+ww); else Mul*=(gnt+ww)/(gnt-ww);
    }
    return Mul;
}

__complex__ double Vainshtein::eSw (int index)
{
    __complex__ double Sw;
    double w0,Xw=0,Yw=0,ww,nu,nu0,OmegaNu0;
    w0=sqrt(pow(Kappa,2)-pow(VMN[index],2));
    nu0=sqrt(pow(Kappa,2)-pow(w0,2));
    OmegaNu0=(index+1)*M_PI*nu0;
    gsl_complex gslBuf;
    int intN=1000;
    // Xw
    for (int i=1; i<intN; i++)
    {
        nu=Kappa*sin((double)i*M_PI/2/intN);
        if (fabs(nu-nu0)<=0.1) Xw+=Omega(nu)/(nu+nu0);
        else Xw+=(Omega(nu)*nu-OmegaNu0)/(pow(nu,2)-pow(nu0,2));
    }
    Xw*=-w0/intN;
    // Yw
    Yw=0;
    double wMax;
    wMax=sqrt(pow(VMN[M_MAX-1],2)-pow(Kappa,2));
    for (int i=1; i<intN; i++)
    {
        ww=(double)i*1000/intN;
        nu=sqrt(pow(Kappa,2)+pow(ww,2));
        Yw+=Omega(nu)/(pow(w0,2)+pow(ww,2));
    }
    Yw*=1000*2/intN;
    Yw+=w0*2/M_PI-(index+1)*M_PI;
    // result
    GSL_SET_COMPLEX(&gslBuf,Xw/2,Yw/2);
    gslBuf=gsl_complex_exp(gslBuf);
    Sw=GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf);
    return Sw; // checked if Yw==Yw_T then ww=Kappa
}

__complex__ double Vainshtein::eSwT (int index)
{
    __complex__ double SwT;
    double w0,Xw=0,Yw=0,ww,nu,nu0,OmegaNu0;
    w0=sqrt(pow(Kappa,2)-pow(MMN[index],2));
    nu0=sqrt(pow(Kappa,2)-pow(w0,2));
    OmegaNu0=index*M_PI*nu0;
    gsl_complex gslBuf;
    int intN=1000;
    // Xw
    for (int i=1; i<intN; i++)
    {
        nu=Kappa*sin((double)i*M_PI/2/intN);
        if (fabs(nu-nu0)<=0.1) Xw+=OmegaT(nu)/(nu+nu0);
        else Xw+=(OmegaT(nu)*nu-OmegaNu0)/(pow(nu,2)-pow(nu0,2));
        //else Xw+=OmegaT(nu)*nu;
    }
    Xw*=-w0/intN;
    Yw=0;
    for (int i=0; i<intN; i++)
    {
        ww=(double)i*1000/intN;
        nu=sqrt(pow(Kappa,2)+pow(ww,2));
        //nu=OmegaT(nu)/(pow(w0,2)+pow(ww,2));
        nu=OmegaT(nu);
        Yw+=nu;
    }
    Yw*=1000*2/intN;
    GSL_SET_COMPLEX(&gslBuf,0,Yw/2);
    Yw+=w0*2/M_PI-index*M_PI;
    // result
    GSL_SET_COMPLEX(&gslBuf,Xw/2,Yw/2);
    //GSL_SET_COMPLEX(&gslBuf,Xw/2,0);
    gslBuf=gsl_complex_exp(gslBuf);
    SwT=GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf);
    return SwT;
}

double Vainshtein::ArgOm(double arg)
{
    gsl_complex gslBuf;
    GSL_SET_COMPLEX(&gslBuf,gsl_sf_bessel_Jn(mI,arg),gsl_sf_bessel_Yn(mI,arg));
    arg=gsl_complex_arg(gslBuf);
    if (arg > M_PI/2) arg-=M_PI;
    if (arg < -M_PI/2) arg+=M_PI;
    arg+=M_PI/2;
    return arg;
}

double Vainshtein::Omega(double arg)
{
    if (arg==0) return arg;
    double Om=0;
    int check=0;
    while (arg>VMN[check] && check<5) check++;
    if (check<5) Om=ArgOm(arg)+M_PI*check;
    else
    {
        double lagrange,points[2][4];
        for (int i=0; i<4; i++)
        {
            points[0][i]=(VMN[check+i]+VMN[check+i+1])/2;
            points[1][i]=ArgOm(points[0][i])+M_PI*(check+i+1);
        }
        for (int i=0; i<4; i++)
        {
            lagrange=1;
            for (int j=0; j<4; j++)
            {
                if (j!=i) lagrange*=(arg-points[0][j])/(points[0][i]-points[0][j]);
            }
            Om+=points[1][i]*lagrange;
        }
    }
    return Om;
}

double Vainshtein::ArgOmT(double arg)
{
    gsl_complex gslBuf;
    GSL_SET_COMPLEX(&gslBuf,besselS(mI,arg),mI/arg*gsl_sf_bessel_Yn(mI,arg)-gsl_sf_bessel_Yn(mI+1,arg));
    arg=gsl_complex_arg(gslBuf);
    if (arg > M_PI/2) arg-=M_PI;
    if (arg < -M_PI/2) arg+=M_PI;
    arg-=M_PI/2;
    return arg;
}

double Vainshtein::OmegaT(double arg)
{
    if (arg==0) return arg;
    double Om=0;
    int check=0;
    while (arg>MMN[check] && check<5) check++;
    if (check<5) Om=ArgOmT(arg)+M_PI*check;
    else
    {
        double lagrange,points[2][4];
        for (int i=0; i<4; i++)
        {
            points[0][i]=(MMN[check+i]+MMN[check+i+1])/2;
            points[1][i]=ArgOmT(points[0][i])+M_PI*(check+i+1);
        }
        for (int i=0; i<4; i++)
        {
            lagrange=1;
            for (int j=0; j<4; j++)
            {
                if (j!=i) lagrange*=(arg-points[0][j])/(points[0][i]-points[0][j]);
            }
            Om+=points[1][i]*lagrange;
        }
    }
    return Om;
}

__complex__ double Vainshtein::eeDelta(double w0)
{
    __complex__ double Sw;
    double Xw=0,Yw=0,ww,nu;
    gsl_complex gslBuf;
    int intN=1000;
    // Xw
    for (int i=1; i<intN; i++)
    {
        ww=(double)i*2*Kappa/intN-Kappa;
        nu=sqrt(pow(Kappa,2)-pow(ww,2));
        Xw+=(Omega(nu)-OmegaT(nu))/(ww-w0);
    }
    Xw*=Kappa*2/intN/M_PI;
    // Yw
    Yw=0;
    intN=1000;
    double wMax,Yw1=0;
    wMax=sqrt(pow(VMN[M_MAX-1],2)-pow(Kappa,2));
    for (int i=1; i<intN; i++)
    {
        ww=(double)i*2*wMax/intN-wMax;
        nu=sqrt(pow(Kappa,2)+pow(ww,2));
        Yw+=Omega(nu)/(pow(w0,2)+pow(ww,2));
    }
    Yw*=wMax*2/intN*w0/M_PI;
    wMax=sqrt(pow(MMN[M_MAX-1],2)-pow(Kappa,2));
    for (int i=1; i<intN; i++)
    {
        ww=(double)i*2*wMax/intN-wMax;
        nu=sqrt(pow(Kappa,2)+pow(ww,2));
        Yw1+=OmegaT(nu)/(pow(w0,2)+pow(ww,2));
    }
    Yw1*=wMax*2/intN*w0/M_PI;
    Yw-=Yw1;
    // result
    GSL_SET_COMPLEX(&gslBuf,Xw/2,Yw/2);
    gslBuf=gsl_complex_exp(gslBuf);
    Sw=GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf);
    return Sw;
}


